subroutine init_react_group(lunit)
  use param
  use group_param,only:react_group
  use group_param,only:n_react_group
  use label_data_param,only:obj_type
  use label_data_param,only:n_obj_type
  use react_data_param,only:react_type
  use react_data_param,only:n_react_type
  use sink_data_param,only:sink_type
  use sink_data_param,only:n_sink_type
  use geom_data_param,only:geom_obj
  use geom_data_param,only:n_geom_obj
  implicit none

  interface
     function reallocInt(p,n)
       integer,pointer,dimension(:) :: reallocInt,p
       integer,intent(in) :: n
     end function reallocInt
  end interface

  integer,intent(in) :: lunit

  integer :: i,j,k
  integer :: n,nr,foo
  integer :: gobj1,gobj2,lid1,lid2,gid1,gid2,id

  if (n_react_group.eq.0) return

  ! obj_type : record the dissociation of the object
  do i=1,n_obj_type
     nr=0
     do j=1,n_react_group
        if (react_group(j)%react_type.eq.2) then
           do k=1,obj_type(i)%n_group_belong
              n=obj_type(i)%group_list(k)
              if (react_group(j)%target1.eq.n) then
                 nr=nr+1
                 obj_type(i)%n_diss_react_group=nr
                 obj_type(i)%diss_react_group=>&
                      reallocInt(obj_type(i)%diss_react_group,nr)
                 obj_type(i)%diss_react_group(nr)=j
                 exit
              end if
           end do
        end if
     end do
  end do

  ! react_type
  do i=1,n_react_type
     nr=0
     do j=1,n_react_group
        if (react_group(j)%react_type.eq.1) then
           foo=0
           lid1=react_type(i)%obj1
           lid2=react_type(i)%obj2
           gid1=react_group(j)%target1
           gid2=react_group(j)%target2
           ! obj1:group1,obj2:group2
           do k=1,obj_type(lid1)%n_group_belong
              n=obj_type(lid1)%group_list(k)
              if (n.eq.gid1) then
                 foo=foo+1
                 exit
              end if
           end do
           do k=1,obj_type(lid2)%n_group_belong
              n=obj_type(lid2)%group_list(k)
              if (n.eq.gid2) then
                 foo=foo+1
                 exit
              end if
           end do
           if (foo.eq.2) then
              nr=nr+1
              react_type(i)%n_react_group=nr
              react_type(i)%react_group=>&
                   reallocInt(react_type(i)%react_group,nr)
              react_type(i)%react_group(nr)=j
              cycle
           end if
           ! obj1:group2,obj2:group1
           foo=0
           do k=1,obj_type(lid1)%n_group_belong
              n=obj_type(lid1)%group_list(k)
              if (n.eq.gid2) then
                 foo=foo+1
                 exit
              end if
           end do
           do k=1,obj_type(lid2)%n_group_belong
              n=obj_type(lid2)%group_list(k)
              if (n.eq.gid1) then
                 foo=foo+1
                 exit
              end if
           end do
           if (foo.eq.2) then
              nr=nr+1
              react_type(i)%n_react_group=nr
              react_type(i)%react_group=>&
                   reallocInt(react_type(i)%react_group,nr)
              react_type(i)%react_group(nr)=j
              cycle
           end if
        end if
     end do
  end do

  ! geom_obj
  do i=1,n_geom_obj
     id=geom_obj(i)%index
     nr=0
     do j=1,n_react_group
        if (react_group(j)%react_type.eq.4) then
           if (react_group(j)%target1.eq.id) then
              nr=nr+1
              geom_obj(i)%n_diss_react_group=nr
              geom_obj(i)%diss_react_group=>&
                   reallocInt(geom_obj(i)%diss_react_group,nr)
              geom_obj(i)%diss_react_group(nr)=j
           end if
        end if
     end do
  end do

  ! sink_type
  do i=1,n_sink_type
     nr=0
     do j=1,n_react_group
        if (react_group(j)%react_type.eq.3) then
           foo=0
           gobj1=react_group(j)%target1
           gobj2=sink_type(i)%gid
           gid1=react_group(j)%target2
           lid1=sink_type(i)%lid
           if (gobj1.ne.gobj2) cycle
           do k=1,obj_type(lid1)%n_group_belong
              n=obj_type(lid1)%group_list(k)
              if (gid1.eq.n) then
                 nr=nr+1
                 sink_type(i)%n_react_group=nr
                 sink_type(i)%react_group=>&
                      reallocInt(sink_type(i)%react_group,nr)
                 sink_type(i)%react_group(nr)=j
                 exit
              end if
           end do
        end if
     end do
  end do

  ! escape thru the surface
  do i=1,n_obj_type
     nr=0
     do j=1,n_react_group
        if (react_group(j)%react_type.eq.5) then
           do k=1,obj_type(i)%n_group_belong
              n=obj_type(i)%group_list(k)
              if (react_group(j)%target1.eq.n) then
                 nr=nr+1
                 obj_type(i)%n_free_react_group=nr
                 obj_type(i)%free_react_group=>&
                      reallocInt(obj_type(i)%free_react_group,nr)
                 obj_type(i)%free_react_group(nr)=j
                 exit
              end if
           end do
        end if
     end do
  end do

  ! Initialize the members of the react_group
  do i=1,n_react_group
     react_group(i)%a_count=0.0
     react_group(i)%a_rate=0.0
  end do

  return

end subroutine init_react_group
